/* Copyright 2012 Dorival de Moraes Pedroso. All rights reserved.
   Use of this source code is governed by a BSD-style
   license that can be found in the LICENSE file. */

#include <iostream>
#include <cstdlib>
#include <cstdio>

#include "mesh/smesh.h"
#include "mesh/umesh.h"
#include "mesh/enrich.h"

extern "C" {

int addlines(GemLab::Mesh *mesh, long ndim,
        long nle,
        long *le_rtags, long   *le_jtags, long *le_quads,
        long le_npts,   double *le_pts,
        long *le_nconn, long   *le_conn,  long *le_btags) {
    if (nle < 1) return 0;
    try {
        for (long i=0; i<nle; ++i) {
            // coordinates
            GemLab::Array<GemLab::Vec_t> verts(le_nconn[i]);
            for (long j=0; j<le_nconn[i]; ++j) {
                long p = le_conn[i*le_nconn[i] + j];
                verts[j].resize(ndim);
                for (long k=0; k<ndim; ++k) {
                    verts[j](k) = le_pts[p*ndim + k];
                }
            }
            // boundary tags
            GemLab::Array<int> btags(2);
            btags[0] = le_btags[i*2 + 0];
            btags[1] = le_btags[i*2 + 1];
            // generate
            GemLab::EnrichMeshBisect(mesh, verts, btags, le_rtags[i], le_jtags[i], (le_quads[i] > 0 ? true : false));
        }
        return 0;
    }
    GEMLAB_CATCH;
    return 1;
}

int umeshgen(char const *fn, long quad, long ndim, long nparts, long pfull,
        double globalmax,
        long *tags,  double *maxs, long nsef, long *sftags,
        long nreg,   double *regs,
        long nhol,   double *hols,
        long npts,   double *pts,
        long *nconn, long   *conn,
        long nle,
        long *le_rtags, long   *le_jtags,  long   *le_quads,
        long le_npts,   double *le_pts,
        long *le_nconn, long   *le_conn,   long   *le_btags,
        long nvt,
        long *vt_tags,  double *vt_coords, double *vt_tols) {
    try {

        // number of regions, holes, points, segments/facets
        GemLab::UMesh mesh(ndim);
        mesh.Set(npts, nsef, nreg, nhol);

        // regions
        for (long i=0; i<nreg; ++i) {
            if (ndim == 2) {
                mesh.SetReg(i, tags[i], maxs[i], regs[i*ndim+0], regs[i*ndim+1]);
            } else {
                mesh.SetReg(i, tags[i], maxs[i], regs[i*ndim+0], regs[i*ndim+1], regs[i*ndim+2]);
            }
        }

        // holes
        for (long i=0; i<nhol; ++i) {
            if (ndim == 2) {
                mesh.SetHol(i, hols[i*ndim+0], hols[i*ndim+1]);
            } else {
                mesh.SetHol(i, hols[i*ndim+0], hols[i*ndim+1], hols[i*ndim+2]);
            }
        }

        // points
        for (long i=0; i<npts; ++i) {
            if (ndim == 2) {
                mesh.SetPnt(i, 0, pts[i*ndim+0], pts[i*ndim+1]);
            } else {
                mesh.SetPnt(i, 0, pts[i*ndim+0], pts[i*ndim+1], pts[i*ndim+2]);
            }
        }

        // segments/facets
        for (long i=0; i<nsef; ++i) {
            if (ndim == 2) {
                mesh.SetSeg(i, sftags[i], conn[i*nconn[i]+0], conn[i*nconn[i]+1]);
            } else {
                GemLab::Array<int> cc(nconn[i]);
                for (long j=0; j<nconn[i]; ++j) {
                    cc[j] = conn[i*nconn[i] + j];
                }
                mesh.SetFac(i, sftags[i], cc);
            }
        }

        // generate
        bool quiet = true;
        mesh.Generate((quad > 0 ? true : false), globalmax, quiet);

        // add lines (if any)
        addlines(&mesh, ndim, nle, le_rtags, le_jtags, le_quads, le_npts, le_pts, le_nconn, le_conn, le_btags);

        // partition domain
        if (nparts > 0) {
            mesh.PartDomain(nparts, pfull);
        }

        // tag vertices (if any)
        if (ndim == 2) {
            for (long i=0; i<nvt; ++i) {
                mesh.TagVertex(vt_tags[i], vt_coords[i*ndim+0], vt_coords[i*ndim+1], 0, vt_tols[i]);
            }
        } else {
            for (long i=0; i<nvt; ++i) {
                mesh.TagVertex(vt_tags[i], vt_coords[i*ndim+0], vt_coords[i*ndim+1], vt_coords[i*ndim+2], vt_tols[i]);
            }
        }

        // write files
        if (ndim == 2) {
            mesh.WritePLY(fn, false);
            mesh.WriteMPY(fn, true, false);
        } else {
            mesh.WritePLY(fn, true);
        }
        mesh.WriteVTU(fn);
        mesh.WriteMSH(fn);
        return 0;
    }
    GEMLAB_CATCH;
    return 1;
}

int smeshgen(char const *fn, long quad, long ndim, long nparts, long pfull,
        long   nreg,   long   *tags,
        long   *nxs,   long   *nys,  long   *nzs,
        double *axs,   double *ays,  double *azs,
        long   *nlx,   long   *nly,  long   *nlz,
        long   npts,   double *pts,
        long   *nconn, long   *conn, long   *btags,
        long   nle,
        long   *le_rtags, long   *le_jtags,  long   *le_quads,
        long   le_npts,   double *le_pts,
        long   *le_nconn, long   *le_conn,   long   *le_btags,
        long   nvt,
        long   *vt_tags,  double *vt_coords, double *vt_tols) {
    try {

        // blocks
        long nbtag = 2 * ndim;
        GemLab::Array<GemLab::Block> blocks(nreg);
        for (long i=0; i<nreg; ++i) {
            // coordinates
            blocks[i].PreInit(ndim, tags[i], nconn[i]);
            for (long j=0; j<nconn[i]; ++j) {
                long p = conn[i*nconn[i] + j];
                for (long k=0; k<ndim; ++k) {
                    blocks[i].C(k,j) = pts[p*ndim + k];
                }
            }
            blocks[i].PostInit(nconn[i]);
            // divisions
            blocks[i].SetNx(nxs[i], axs[i], (nlx[i] > 0 ? true : false));
            blocks[i].SetNy(nys[i], ays[i], (nly[i] > 0 ? true : false));
            blocks[i].SetNz(nzs[i], azs[i], (nlz[i] > 0 ? true : false));
            // boundary tags
            for (long j=0; j<nbtag; ++j) {
                blocks[i].BryTags[j] = btags[i*nbtag + j];
            }
        }

        // generate
        GemLab::SMesh mesh(ndim);
        mesh.Generate(blocks, (quad > 0 ? true : false));

        // add lines (if any)
        addlines(&mesh, ndim, nle, le_rtags, le_jtags, le_quads, le_npts, le_pts, le_nconn, le_conn, le_btags);

        // partition domain
        if (nparts > 0) {
            mesh.PartDomain(nparts, pfull);
        }

        // tag vertices (if any)
        if (ndim == 2) {
            for (long i=0; i<nvt; ++i) {
                mesh.TagVertex(vt_tags[i], vt_coords[i*ndim+0], vt_coords[i*ndim+1], 0, vt_tols[i]);
            }
        } else {
            for (long i=0; i<nvt; ++i) {
                mesh.TagVertex(vt_tags[i], vt_coords[i*ndim+0], vt_coords[i*ndim+1], vt_coords[i*ndim+2], vt_tols[i]);
            }
        }

        // write files
        if (ndim == 2) {
            mesh.WriteMPY(fn, true, false);
        }
        mesh.WriteVTU(fn);
        mesh.WriteMSH(fn);
        return 0;
    }
    GEMLAB_CATCH;
    return 1;
}

} // extern "C"
